Fluctuation induces evolutionary branching in a modeled microbial ecosystem 
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The impact of environmental fluctuation on species diversity is studied with a model of the evo- 
lutionary ecology of microorganisms. We show that environmental fluctuation induces evolutionary 
branching and assures the consequential coexistence of multiple species. Pairwise invasibility anal- 
ysis is applied to illustrate the speciation process. We also discuss how fluctuation affects species 
diversity. 
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Fluctuation is ubiquitous in nature. Biological sys- 
tems are always exposed to environmental fluctuations; 
hence, the systems have evolved under these fluctuations 
Ql Q- E[ • Ecological systems are no exception 0, [H[ . Di- 
versity, which is one of the most essential properties in 
ecology, has also evolved in the presence of fluctuations. 
The association between diversity and fluctuation has 
been discussed, and concepts such as intermediate dis- 
turbance hypothesis p have been proposed. However, 
there are few discussions about the mechanisms by which 
fluctuation is associated with diversity. 

In this paper, we study a model of the evolutionary 
ecology of microorganisms in the presence of a tempo- 
ral fluctuation. We consider competition for a single re- 
source, one of the simplest but most thoroughly investi- 
gated ecological situations @, • In this case it is often 
believed that only one species can survive by following 
the competitive exclusion principle [9[. However, here 
we demonstrate that temporal fluctuation in the envi- 
ronment induces evolutionary branching and the stable 
coexistence of multiple species. The results and analy- 
sis of the evolutionary dynamics reveal how fluctuation 
facilitates species diversity. 

We adopted a microbial ecosystem because experimen- 
tally based quantitative descriptions of its growth kinet- 
ics are available [Io|, EH • The growth rate limited by a 
resource is described by a saturation function called the 
Monod function [To| . 



/x(s;/x,fe) = jl- 



(1) 



jl oc log(&). The former denotes a trade-off relation: fast 
growth in rich media (large jl) vs. the ability to grow in 
a wide range of resource concentrations (small k). Fol- 
lowing their work, we introduced a continuous genetic 
parameter g in our model which specifies the parameters 
of the Monod function. We chose 
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with an appropriate scale conversion. The mutation in- 
dicates the small change in g. At the population level it 
is described with the diffusion process in g space. 

Now we give the model equation by referring to the 
chemostat experiments [7j. The individual genotype 
grows at its own growth rate /i(s; jl g , k g ) and decreases 
at the same rate 7. The resource s is supplied with the 
time-dependent function c(t) and consumed by all popu- 
lations in proportion to their growth rates. Let x g be the 
population density of a genotype. The model is given by 
a partial differential equation with globally interacting 
variable s, 
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We choose D = 5 • 10 -7 and 7 = 0.11 in the following 
simulations. 

As the feeding manner, we choose a continuous supply 
and additions with periodic pulses 



c(t) = Co + Ci 



(4) 



where s denotes the concentration of the limited resource. 
The two parameters, the maximum growth rate jl and 
the half-saturation constant k, are genetically determined 
and characterize the strategy for the resource utilization 
of the genotype. 

Moreover, recent experimental developments allow us 
to observe the evolution and ecological diversity of mi- 
croorganisms directly H, [TH, [l3j]. The data in [ll[ in- 
dicate the relations between parameters of the Monod 
function among different genotypes, such as the posi- 
tive correlation djl/dk > 0, and the logarithmic relation 



The periodic additions represent the input fluctuation. In 
this paper we choose T = 5. We checked that the results 
are reproduced over a wide range of the period. Another 
parameter, A, is introduced to describe the intensity of 
the fluctuation. The amount of supplied resource is given 
as Co = 1 — A and c\ = A • T, where the total amount of 
supply in a period is fixed. 

Since we introduced the continuous genotype space, 
the population is distributed in the space, and the geno- 
type and species do not hold a one-to-one correspon- 
dence. We use the term quasi-species [14[ for a group 
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FIG. 1: Evolutionary processes as the pattern dynamics of 
population densities in the genotype spaces. The conditions 
of three different intensities of fluctuation A = 0.5, 0.7, and 
0.9 are shown. The vertical axis denote the genotype space, 
and the population density normalized in time with period 
T = 5 is plotted with gray-scale in a logarithmic scale. Ini- 
tially, static conditions (co, ci) = (1— A, 0.0) are set, and stable 
genotype compositions are provided. Fluctuations are intro- 
duced at t = 5000 (co,ci) = (1 — A, 5 • A), and evolutionary 
branching is clearly observed after that. 

of the population which has a unimodal distribution in 
the genotype space, and is represented by the genotype 
at the peak of the distribution. If the mutation rate is 
sufficiently small, the population can be approximately 
replaced by the population of the representative geno- 
type. 

In the simulations, population densities less than 10 -8 
are replaced with zero, which denotes the discreteness of 
the population. 

Fig. [1] shows the evolutionary dynamics in the three 
different conditions, A = 0.5, 0.7, and 0.9. In each simula- 
tion we initially supply the resource only in a continuous 
way and the genotype composition of the resident popu- 
lation shows a unimodal distribution in g space; i.e., one 
quasi-species exists. This agrees with the competitive 
exclusion principle. However, introducing fluctuation in 
the resource supply (the timings are marked with arrows 
in Fig. [1]), we find dynamic changes of the genotype dis- 
tributions which lead to evolutionary branching. Each of 
the emerging branches has a unimodal shape and has no 
connection with the others. Therefore, each of them is 
regarded as an independent quasi-species. This implies 
that the coexistence of multiple quasi-species is attained. 

When the periodic supply is introduced, the genotype 
distribution changes through two phases, gradual evolu- 
tion and branching (indicated by I and II in Fig. [l]-a). In 
the first phase, the distribution moves gradually to higher 
g in the genotype space, keeping a unimodal shape, which 
indicates the gradual evolution of the quasi-species. The 
branching starts after the movement stops. 

The final states of the three figures in Fig. [I] have sim- 
ilar genotype compositions. Therefore, the states with 
two quasi-species are robust against the intensity of fluc- 
tuation A. However, features of evolutionary branching 




FIG. 2: (color online) The dependences of characteristic 
genotypes on the intensity of fluctuation A. Data from the 
evolutionary simulations are plotted with lines. The represen- 
tative genotypes at the end of the first phase are plotted with 
dashed lines, and the representative genotypes at the final 
state are plotted with continuous lines. Characteristic values 
of PIPs (discussed below), g c and g ex , are plotted with red 
squares and green circles, respectively. g c values show good 
correspondences with genotypes at the end of first phase, g e x 
values are related with genotypes at coexisting states. 



dynamics are different. The differences are mainly char- 
acterized by the states at the onset of branching (the 
states at the end of the first phase). In Fig. [2j we plot 
lines of the genotype representing the quasi-species at the 
end of the first phase and those of the genotype (s) rep- 
resenting the quasi-species at the final states against A. 
Branching was observed in the range A ~ 0.40 ~ 0.96. At 
the no-branching conditions, only the first phase appears, 
and the two lines coincide. 

The oscillatory time series of the system after it reaches 
the final stationary state for A = 0.9, and the Monod 
function for the resident quasi-species are shown in Fig. 
[3j At this state, the two quasi-species have different 
strategies based on the trade-off relation. The difference 
in the available ranges of resource concentration repre- 
sents the difference of strategies: a strategy adapted for 
high concentration and a strategy adapted for a wide 
range of concentration. Therefore, their coexistence can 
be regarded as a separation of niches in the resource con- 
centration space. The temporal fluctuation enables them 
to utilize the full extent of the concentration space. 

In the following, we investigate the mechanism of evo- 
lutionary branching dynamics and discuss how fluctua- 
tion changes the utility of environmental states. 

We start the analysis with the nonfluctuating condi- 
tion, which corresponds to the situations before the pe- 
riodic fluctuations are added in the three simulations in 
Fig. [TJ In this condition, one quasi-species is resident. 

Suppose a system containing one genotype g = a with- 
out mutation. The resource concentration comes to an 
equilibrium s a which holds fi(s a ; fta,k a ) — 7 = 0. Note 
that s g depends only on the traits of the genotype and 
is independent from the amount of supply cq. Thus s g 
is the value characteristic to the genotype. If another 
genotype g = /3 with a smaller equilibrium (sp < s a ) 
is introduced into the state, the population of f3 grows 
(/i(s a ; ftp, kp)—j > 0). This raises the total consumption 
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FIG. 3: (color online) (a) Oscillation of the resource con- 
centration (solid line) and populations of the quasi-species 
(dashed lines) for A = 0.9. Data are taken after the branch- 
ing dynamics are completed, (b) The Monod function for 
genotypes representing these quasi-species. The fluctuation 
range of resource concentration is indicated. In this range, the 
faster- growing species change as a consequence of the trade-off 
relation. 



of the resource, and the resource concentration decreases. 
Then the genotype a decreases and the system comes to 
a new equilibrium state with the population of genotype 
ft and the resource concentration sp. The monotonic in- 
crease of the Monod function guarantees these processes. 

Therefore, the relation between genotypes is deter- 
mined by the values of s g , and the genotype g c which 
has the minimum value is the fittest genotype. The res- 
ident genotypes are replaced with genotypes with lower 
equilibria sequentially until genotype g c is reached, and 
the resource concentration comes to s 9c . In our model 
g c 0.24 corresponds to the representative genotype of 
quasi-species at the initial static conditions in simula- 
tions. 

When fluctuation in the resource supply is introduced, 
the above simple discussion is not directly applicable. 
This is because the growth or decay of populations de- 
pends on oscillation patterns (as an example of the oscil- 
lation profile, see Fig. [3]-a). Such a situation is in con- 
trast to the above condition, where only the static value 
of s determines the growth rate. To analyze the fluctua- 
tion conditions we introduce the pairwise-invasibility plot 
(PIP), which gives the relation between genotypes and 
enables us to describe the evolutionary dynamics with 
the process of invasion and annihilation. The discussions 
are based on Geritz et al. [HI- 

PIP is constructed with the following procedure. Sup- 
pose a system with one genotype g = a, which is given 

by 



s = c(t) — /i(s; [la-, ka)x a - 



(5) 



Integration gives a periodic oscillation function of re- 
source concentration s a (t) with period T after some re- 
laxation time. If a small population of another genotype 
g = ft is introduced into the system, the average growth 
rate of ft is calculated by 



a (a, ft) 



1 f T 

T Jo { fJ, Mt)lp>/3,kp)-'y}dt. 
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The sign of cr(a, ft) determines the invasibility of ft into 
the a population. Calculating a (a, ft) for all pairs of 
genotypes and shading the areas of positive a (a, ft) on 
the a- ft plane, PIP is obtained. The PIPs corresponding 
to Fig. [U are shown in Fig. [U 

Here we investigate the structures of the PIPs and use 
them to illustrate evolutionary dynamics. The PIPs are 
constructed with two boundary lines: one is a diagonal 
line a = ft [17], and the other is a curved line 



0\(a,/3) =0, 



(7) 



the shape of which determines the characteristics of PIP. 
First, we refer to the genotype at the intersection be- 
tween lines as g c [18], which is important because of the 
following singular property. Any genotype lower than g c 
is invaded by genotypes higher than it (area above the 
diagonal line is shaded) and any genotype higher than 
g c is invaded by genotypes lower than it (area below the 
diagonal line is shaded). This indicates that, as part of 
the evolutionary process, a resident genotype is replaced 
by invading genotypes closer to g c one after another, and 
the genotype converges to g r ,. Therefore g c is called the 
convergent stable genotype [l5[. The convergent process 
corresponds to the gradual evolution of the quasi-species 
at the first phase in simulations. It stops when the rep- 
resentative genotype agrees with g c . The agreement be- 
tween them is shown in Fig. O 

In Fig. [4j the vertical line on g c passes through the 
shaded regions. This means that there are genotypes 
invasible for the population of g Cl and this invasibility 
promotes evolutionary branching. At the end of the first 
phase the population distributes around g c in the geno- 
type space. If the tail of the distribution covers the region 
of invasible genotypes, the population of these genotypes 
grow and form another branch. Therefore, evolutionary 
branching succeeds as the second phase. 

The features of branching dynamics are characterized 
by the relation between g c and genotypes invasible to g c . 
In PIP for A = 0.5 genotypes higher than g c are invasible 
to g Cl both higher and lower genotypes are invasible to 
g c for A = 0.7, and lower genotypes are invasible to g c for 
A = 0.9. Correspondingly in Fig. [TJ the upper branch 
is born from the lower one for A = 0.5, the unimodal 
shape separates symmetrically for A = 0.7, and the lower 
branch is born from the higher one for A = 0.9. 

In addition, the resident quasi-species at the final state 
are related with genotypes g ex , which give the local max- 
imum or minimum for a in eq. Q, 



d/3, 
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0. 
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These genotypes can be the locally fittest genotype in the 
sense that they can invade into the maximum ranges of 
genotypes (shaded regions bounded by eq. (0) have local 
maximum widths at the points). We plot g ex in Fig. [2j In 
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FIG. 4: (color online) Pairwise invasibility plot for A = 
0.5, 0.7, and 0.9 (corresponding to Fig. [TJ. g c and g ex are plot- 
ted with red squares and green squares, respectively. Vertical 
lines on g c (red solid lines) pass through the shaded regions, 
which denotes genotypes invasible to g c . 



the branching range of A, the representative genotypes for 
the final states stay in the vicinities of g ex . In particular, 
at the onset of the branching range, one quasi-species 
coincides with g c and the other coincides with g ex . This 
is because at the onset, the curved boundary line 6\ is 



tangent to the vertical line on g c at (3 



and an 



additional quasi-species arises at the tangent point. 

In summary, evolutionary branching dynamics induced 
by environmental fluctuation is reported in a model of a 
microbial ecosystem competing for a single resource. A 



pairwise-invasibility plot, which was introduced in [15] 
for a static environment with multiple parameters, is ex- 
tended to a fluctuating environment and is applied here 
to illustrate evolutionary dynamics. 

Previous studies Q have reported the coexistence 
of two species with given parameters in the presence of 
fluctuation. However, the occurrence of the branching 
and evolutionary stability of the coexisting state have 
remained open questions. Here, we give clear demon- 
strations of them. 

Although we use a model of a microbial ecosystem, the 
branching mechanism we found does not depend on the 
detail of the model. Therefore, the phenomena must be 
found in more general classes of ecosystems. 

MacArthur and Levins [l6[ indicated that the num- 
ber of independently adjustable environmental parame- 
ters corresponds to the maximum number of coexisting 
species as an extended and refined version of the com- 
petitive exclusion principle. Applying it to our system 
for the static case, resource concentration is only the ad- 
justable parameter, and it leads to the existence of only 
one quasi-species. However, the stable co-existence of 
two quasi-species is seen in cases where fluctuation is 
present. This suggests that the number of adjustable pa- 
rameters is increased by introducing environmental fluc- 
tuation. When the system is always in an oscillatory 



state, whether one can grow or not depends on the com- 
prehensive details of the oscillation profile of the param- 
eter. In other words, the whole oscillation profile is the 
environmental state to be adjusted. The variety of the 
profile can be described with several dimensional param- 
eters which forms a subspace of a functional space. In 
this way, introducing fluctuation expands the dimensions 
of adjustable environmental parameters, and it enables 
branching and the coexistence of genotypes with differ- 
ent strategies. The diversity of oscillation profiles is not 
restricted to two-dimensional parameter space. There- 
fore, by the principle, the coexistence of more than two 
species is possible. Although oscillation is introduced as 
an input fluctuation in this work, the self-sustained oscil- 
lation in population dynamics, if one supposes a system 
of generating it, can also induce evolutionary branching 
based on the same mechanism. 

The author is grateful to K. Kaneko, S. Ishihara, K. 
Fujimoto, and M. Inoue for their helpful suggestions. 
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